# vizualize spatial sst dynamics 
library("raster")
library("tidyverse")
library("rasterVis")

sst <- stack("input/cmems_mod_blk_phy-tem_my_2.5km_P1M-m_1715770301759.nc")


## calc year anomaly
sst19931995 <- sst[[1:36]]
sst20192021 <- sst[[313:349]]

avg1995 <- mean(sst19931995)
avg2021 <- mean(sst20192021)

# sst avg 1993-1995 vs 2019-2021 
at <- c(10, 12, 14, 16, 18, 20)
cc <- list(at=at, ## where the colors change
           labels=list(
             at=at## where to print labels
           ))
png("output/mapsst/annual-sst-1993-1995.png", res=180, width=1200, height=900)
print({
  levelplot(avg1995,  contour = F, region = T, margin = F, main = paste0("Среднегодовая ТПВ 1993-1995"), par.settings = magmaTheme, at=at, colorkey = cc)
})
dev.off()

png("output/mapsst/annual-sst-2019-2021.png", res=180, width=1200, height=900)
print({
  levelplot(avg2021,  contour = F, region = T, margin = F, main = paste0("Среднегодовая ТПВ 2019-2021"), par.settings = magmaTheme, at=at, colorkey = cc)
})
dev.off()

# annual year anomaly
# sst avg 1993-1995 vs 2019-2021 
anomaly <- avg2021 - avg1995
at <- c(0, 0.5, 1, 1.5, 2.0, 2.5, 3.0, 3.5)
cc <- list(at=at, ## where the colors change
           labels=list(
             at=at## where to print labels
           ))
png("output/mapsst/anomaly-sst-anual-1993-to-2021.png", res=180, width=1200, height=900)
print({
  levelplot(anomaly,  contour = F, region = T, margin = F, main = paste0("Аномалии ТПВ 2019-2021 относительно 1993-1995"), par.settings = BuRdTheme, at=at, colorkey = cc)
})
dev.off()

## calc winter anomaly
sstwint19941996 <- sst[[c(12,13,14,24,25,26,36,37,38)]]
sstwint20192021 <- sst[[c(312,313,314,324,325,326,336,337,338)]]

avgwint19931995 <- mean(sstwint19941996)
avgwint20192021 <- mean(sstwint20192021)

# sst avg 1993-1995 vs 2019-2021 
at <- c(2, 4, 6, 8, 10, 12, 14)
cc <- list(at=at, ## where the colors change
           labels=list(
             at=at## where to print labels
           ))
png("output/mapsst/winter-sst-1994-1996.png", res=180, width=1200, height=900)
print({
  levelplot(avgwint19931995,  contour = F, region = T, margin = F, main = paste0("Среднезимние ТПВ 1994-1996"), par.settings = magmaTheme, at=at, colorkey = cc)
})
dev.off()

png("output/mapsst/winter-sst-2019-2021.png", res=180, width=1200, height=900)
print({
  levelplot(avgwint20192021,  contour = F, region = T, margin = F, main = paste0("Среднезимние ТПВ 2019-2021"), par.settings = magmaTheme, at=at, colorkey = cc)
})
dev.off()

anomalywint <- avgwint20192021 - avgwint19931995
at <- c(0, 0.5, 1, 1.5, 2.0, 2.5, 3.0, 3.5)
cc <- list(at=at, ## where the colors change
           labels=list(
             at=at## where to print labels
           ))
png("output/mapsst/anomaly-sst-winter-1994-to-2021.png", res=180, width=1200, height=900)
print({
  levelplot(anomalywint,  contour = F, region = T, margin = F, main = paste0("Аномалии среднезимних ТПВ 2019-2021 относительно 1993-1995"), par.settings = BuRdTheme, at=at, colorkey = cc)
})
dev.off()

## calc statistics 
print(paste0("Annual anomaly mean=",cellStats(anomaly, stat='mean'), ', sd=', cellStats(anomaly, stat='sd')))
print(paste0("Winter anomaly mean=",cellStats(anomalywint, stat='mean'), ', sd=', cellStats(anomalywint, stat='sd')))

